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Abstract 

Assuming the existence of a Bose-Einstein condensate composed of the majority 
of a sample of ultracold, trapped atoms, perturbative treatments to incorporate the 
non-condensate fraction are common. Here we describe how this may be carried 
out in an explicitly number-conserving fashion, providing a common framework 
for the work of various authors; we also briefly consider issues of implementation, 
validity and application of such methods. 

1 Introduction 

Since the first successful experiments in observing Bose-Einstein condensation (BEC) 
in cold trapped atoms [ 1 , 2|, there have been dramatic advances into strongly interact- 
ing (via Feshbach resonances Q) and strongly correlated regimes (e.g., within optical 
lattices f?, '5l|). Nevertheless, one can still speak of a 'typical' BEC experiment as 
consisting of a weakly interacting gas of alkali atoms, held within a confining poten- 
tial formed by laser or magnetic fields, where the condensate fraction incorporates the 
substantial majority of constituent atoms. In this chapter we will describe perturbative 
approaches, based around the existence of such a significant condensate fraction, in 
such a way that the many-body system is in a number eigenstate (number conserving) 
rather than a coherent state (symmetry breaking) The material presented here is in- 
tended to be a systematized amalgam of the presentations of C. W. Gardiner |6 1, Castin 
and Dum Q, and of S. A. Gardiner and Morgan |8 1, with some additional observations. 

We begin by contrasting so-called number-conserving approaches with the more 
conventional assumption of symmetry-breaking, also addressing the motivation for 
considering a number-conserving alternative, before describing number-conserving equiv- 
alents to the quadratic Bogolibov Hamiltonian ||9][T0]|6][lT| and the linearized Bogoliubov- 
de Gennes equations jT] [121 |T3l . We then cover extending this approach to a second- 
order minimal self-consistent treatment of dynamics, before concluding with consider- 
ations of implementation, validity and application. Throughout, due to the preponder- 
ance of explicit time-dependences, t arguments will only appear when a new quantity 
is introduced, as appropriate. 
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2 Methodology: A Number-Conserving Perturbative Ap- 
proach 



2.1 Number- Conserving versus Symmetry-Breaking Approaches 

Our starting point is the bosonic binary interaction Hamiltonian 

J t/r4''(r)^o(r)4'(r) + ^ J c/r^^'Xr)*' (r)4'(r)4'(r), (1) 

where the field operators obey [4'(r),4'^(r')] - S(r - r'), and ^o(r) = -{fp- llnijV^ + 
y(r); y(r) is an external potential, m is the atomic mass, and g - A-nJfa/m, with a 
the s-wave scattering length. The binary interaction here is characterised by a contact 
term g6(r - r'), with the conditions for which a renormalization of the consequent 
ultraviolet divergences is necessary described elsewhere fTT' TSl. The number operator 
N - J (ir*l''(r)4'(r) commutes with H; particle number is therefore conserved, and 
stationary states of H must also be eigenstates of .^V. 

If (as is commonly assumed) one defines the condensate as 0(r, f) - (*J'(r)) (the 
mean field) [161 . the field operator can be written 

4'(r) = ,^(r,f) + 'A'(r,f), (2) 

with fluctuations defined by i^'(r) = *J'(r) - (*J'(r)). A scalar function, 0(r) has ev- 
erywhere a well-defined phase, and the U(l) gauge symmetry of the system is bro- 
ken. Consequently, for (*l'(r)) ?t 0, the system must be in a coherent state, i.e., a 
coherent superposition of [TT','?!. The appropriate statistical ensemble is then grand- 
canonical rather than canonical, and one should work with the grand-canonical Hamil- 
tonian K - H - jiN 1 17, 6|. Assuming small fluctuations, approximate expressions 
may then be obtained perturbatively around purely mean-field results. 

Requiring the state of the system to be in a number eigenstate implies (*l'(r)) = 0. 
Data is often statistically averaged from repeated experimental runs, which are unlikely 
to have identical shot-to-shot particle numbers — nevertheless, if each run is for a 
definite (even if unknown) particle number, (*I'(r)) is still = (one has an incoherent 
statistical ensemble of number eigenstates, rather than a coherent superposition). The 
alternative taken 17] |6l |9] [lOl [14l [81 is to select a condensate mode <p^{r, t): 

4'(r) = flo(f)0''(r,f) + <A''(r,f), (3) 

where tff'^(r) is explicitly orthogonal to ^'^(r); formally tff'^ir, t) - J dr'Q(r, r', f)4'(r'), 
with the projector defined by (3(r,r',f) - 6(r-r')-(f)^(r,t)(f)^*{r' ,t). Following Penrose 
and Onsager ITSl . one can use the single-body density matrix p(r, r', f) - (4*' (r')*P(r)) 
to formally define the condensate mode as the instantaneous eigenfunction of p(r, r', t) 
with the largest eigenvalue Ndt) (the condensate number): 

J dr'p(r, r', t)^^(r', t) = A^c(0<^^(r, 0- (4) 

It follows directly that Nc - (a^^ao), and (fly0''^(r)) = 17|. Such a prescription is 
consistent with off-diagonal long-range order, i.e., p(r, r', f) = Nc{t)(l>^{r, t)<p^*(r', t) + 
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{iff^'^ir' ,t)^^{r,t)) |fT9l . To maintain mutual orthogonality, it is expected that nonlo- 
cal terms will arise in the dynamical equations, as indeed occurs liZiiSj. The generic 
appearance of nonlocal terms is associated with the orthogonality of condensate and 
non-condensate, however, not number conservation per se. Similarly, although the 
commutation relations [(^'(r), i^'' (r')] - 5(r - r') are exactly bosonic, if the fluctua- 
tion term is defined orthogonal to the condensate mode, its commutation relations are 
projective: [!if^{r),!if^Hr')] =(3(r,r') QEI. 



2.2 Key Concept: A Number-Conserving Fluctuation Operator 

This Chapter is concerned with describing a perturbative approach, closely analogous 
to treatments described in the three preceeding chapters, but leading to a theory which 
is guaranteed to be number-conserving by construction and does not therefore need to 
treat the condensate part as a mean field (see also the stochastic approaches of Part II. 
2). This is achieved by introducing an appropriate operator within a number-conserving 
context, to analogously describe equivalent small fluctuations to those about the mean 
field of a symmetry-breaking treatment. 

Recalhng [Eq. (j4])] that (agi^'^(r)) = 0, we may use a^ifi^ir) as the core of an 
appropriate fluctuation operator; in particular, {a[^\j/^{Y)) is not trivially - 0, as a^^ip^ir) 
still conserves total particle number. Assuming an almost fully condensed system (i.e., 
Ac ~ N), we may rescale this fluctuation term to scale approximately as ^^{v), defining 

A(r,f)= ^^^(OiAV,?), (5) 

such that <A(r)> = and [A(r), A"'"(r')] ~ [^r^ir), ^r^Hr')] = <3(r, r'). Castin and Dum 
define A(r, f) - N^^^^aj^(t)ifr^(r,ty, if the system is in a number eigenstate, and all 
other operators are arranged in number-conserving pairs in any approximate Hamil- 
tonian or equation of motion, .^V behaves exactly as a number A^, as made explicit 
here. The comparable ;f-(r, t) operators of C. W. Gardiner [SI are defined through 
X(r,t) = N-^^^al(t)4^^(r,t), where N = A'Ut)A(t) and A(f) = y/N/NMaoit); hence 
(x(r)} is only approximately zero, and [x(r),x''(r')] - (Nc /N)Q(r, r'). In their approx- 



imate application, to Bogoliubov order (see Section 2.3 below), x(r) and A(r) yield 
identical results. 

An equivalent operator compatible with the treatment of Girardeau and Arnowitt 
ll9l[ip| is Ac(r, t) - Ncity^'^al{t)i//'^(r, f), whereas Gardiner and Morgan use A(r, t) - 
Nc(ty^^^a^Q(t)iJr^(r,t) |8|. The operators Ac(r) and A(r) are scaled by the condensate 
number, and hence intended to be better-tailored to describing larger non-condensate 
fractions. Noting the exact identities [Ac(r), Ac(r')] = <3(r,r') and Ai(r)Ac(r) = 
1^^' (r)0''^(r), A(;(r) may seem an attractive choice; (Ac(r)) is not identically zero, how- 
ever, whereas {A(r)) = is an exact identity. Hence, in the later sections of this chapter 
we use A(r). 



2.3 Number- Conserving Bogoliubov Treatment 

To determine an appropriate Bogoliubov Hamiltonian, while explicitly maintaining 
number-conservation, we first substitute Eq. ([3]) into Eq. ([TJ, neglecting all terms of 
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greater than quadratic order in ^''^(r), i^^^{r). We set gN - U, and express the trun- 
cated Hamiltonian in terms of A(r), A' (r). The operators Nc NdNc - l)/N appearing 
in the naively zeroth-order term in this truncated Hamihonian must be replaced by 
- / dr'4>^\r')^^(r') ~ N - J t/r'A"'"(r')A(r') and A? - 1 - 2 / dr'{(f^Hr')4>^{r') ^ 
N - 1 - 2 Jt/r'A' (r')A(r'), respectively, as these include second-order terms; we also 

set Nc/ "V^ — > Va^, N/Nc — > 1, and smaller terms — > 0. The consistent second-order 
Hamiltonian is thus ||6l 

H2 =//o + Va^ J dr {0''*(r) [^o(r) + Ultp^'irf] A(r) + H.c.} 

+ J drA\r)[ho(r) + 2U\4>''(r)f]Mr) 

+ IJ ^r[0^*(r)2A(r)2+H.c.]-|J rfr|A 

- j dr'AHr')A(r') j flfr0^*(r) \ho(r) + U\cf>^irf] cf>^(r). 

Minimizing Hq/N = j dr(p^* ir)[hn{r) + iU/2)\(l)^ir)\^]<p^{r) with respect to ^^(r) (sta- 
tionary, and constrained to preserve unit norm) yields the time-independent Gross- 
Pitaevskii equation (GPE) 



(6) 



[/2o(r) + f/|<^^(r)|2]^^(r) = ^0^(r), 



(7) 



where A arises as a Lagrange multiplier, and takes the form of a nonlinear eigenvalue 
16]. Substituting Eq. (Ml back into H2 eliminates all terms linear in A(r), A' (r), leaving 
H2 in a quadratic form suitable for diagonalization by Bogoliubov transformation. 

In an explicitly dynamical treatment, we use H2 and [A(r), A' (r')] ~ (3(r, r') ||6j 
El El to generate the equation of motion ihdA(r)/dt - [A(r), i/2] + ifidK{r)ldt. This 
yields WM 



in-k{r)^^fN f dr'Q(r,r' 
at J 



/zo(r') + t/|0^(r')r - in 



dt 



J dr'Qir, r') [/lo(r') + 2f/|0N(r')P] A(r') 



+ U 
-Mr) 



fdr' 



(3(r,r')A"'"(r')0^(r') -<^(r) 



ih- 



df 



^*(r') 



(8) 



A(r') 



J fl'r>^*(r') 



hair') + f/|0"(r' 



N/„'m2 



ifi 



dt 



Ar'). 



Using (A(r)) = yields the time-dependent GPE 

d 



n-^'^ir) = [^„(r) + f/|0^(r)|2 - a\ 0^(r), 



(9) 



where A = /t/r(^'^*(r) [^o(r) + t/|<;*^(r)p - indldt\ 4>^{r) corresponds to an arbitrary 
global phase for (p^{r) Substituting this result back into Eq. (jsjl simplifies the 



'Note that we regain the time-independent GPE for itidip^(r)/dt = 0. 
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expression, which can be combined with its Hermitian conjugate to form QIHl 



A(r) 



dt\ A"'"(r) 
X(r,r') 



-C(r,r') 



A(r') 



AT(r') 

L(r,r') M(r,r') 
M*(r,r') -L*(r,r') 

L(r, r') =5(r - r')[^o(r') + U\(f>'^(r')f- - A] (10) 
+ U J dr'Q(r,r')\cf>ir'ffQ(r',r'), 

M(r,r')=C/ J" <ir'<3(r,r')<^'^(r')^(3*(r',r'), 

i.e., the standard Bogoliubov-de Gennes equations (EDGE) modified only by the ap- 
pearance of nonlocal projector terms. Note that the treatment of C. W. Gardiner IS) 
avoids employing (A(r)) - 0, instead transforming to a 'condensate picture.' Use 
of this identity is well-motivated, however, and yields equivalent results in a simpler 
fashion. 

The X.(r, r') can now be diagonalized |l2l[Hl: 



X(r,r') = J]ei 



vk(r) 



(«*(r'),-v:(r'))- 



^^(r) 



(-y,(r'),«i(r')) 



(11) 



(k is an arbitrary index) where ((p^(r), 0) or (0, <p^*(r)) are also eigenstates, with eigen- 
value = 0. The spectrum is identical to that associated with the EDGE, and the 
eigenstates differ only in that they are explicitly orthogonal to the condensate mode 
I?]. We decompose (A(r), A"'"(r)) as 



A(r) 
At(r) 



k=\ 



Uk(x) 
Vk(r) 



+ b 



«I,(r) 



(12) 



the orthonormality relations J dr[ul,(r)uic(r) - v'^,,(r)vi:(r)] - dkt, J dr[uk'(r)vic(r) — 
Vk'{v)uii{r)] - 0, apply, meaning that the quasiparticle annihilation operators b/i = 
J c/r[M^(r)A(r) - v^(r)A'(r)] with their Hermitian conjugates form a bosonic alge- 
bra, i.e., [bk,bl] = 5kk' (as [A(r), A"'"(r')] ~ <3(r,r') to current order) lEHS). The 
system dynamics are to this order described by the modified EDGE coupled to the 
time-dependent GPE. This can lead to a situation of rapid non-condensate growth, 
but condensate depletion is not accounted for Il20l |2TI . which effectively describes a 
zero-temperature, infinite particle limit fSl; a treatment to the same order using N^^^^ 
(rather than N^^^^) as an asymptotic expansion parameter [14, 8| is therefore func- 
tionally equivalent, as A^c must x N. Assuming a stationary configuration (such that 
first-order terms are eliminated), H2 can now be written in diagonal form: 



H2^Ho-^J dr\4>''(r)f + |j blh - J drMr)f 



(13) 



i.e., we have carried out an equivalent, number-conserving Eogoliubov transformation. 
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2.4 Second- Order Self- Consistent Treatment 



To account for significant thermal or dynamical depletion, we must go to higher or- 
der. Using second-order perturbation theory]^ Morgan determined a fourth-order ap- 
proximate Hamiltonian, for a system at thermal equilibrium |14|, which had a gapless 
excitation spectrum, in accordance with the Hugenholtz-Pines theorem |22|. Here we 
briefly describe the dynamical second-order treatment of S. A. Gardiner and Morgan 
(HI, in terms of A(r); this is the minimal order necessary for the consistent treatment of 
particle transfer between condensate and non-condensate. 

The first term in iMA(r)/dt = iiMNJdt) A(r)/2Nc + N;^'^(iM[aliJr^(r)]/dt) is of 
cubic order, and is dropped |8 1. One can substitute in the first 3 terms of Eq. (A4) in 
Q, and express the result in terms of A(r, f), A"'"(r, f), = gNc, setting Nc/Nc — > 1, 
■\/Nc/Nc — > 1 / '\/Nc and smaller terms to zero. Taking the expectation value then yields 
the generalized GPE (GGPE) 



d 



<A(r)2> U 



J dr'lcp^'ir't 



X f<A^(r')A(r))^*^(r') + <^N*(r')<A(r')A(r))] ., 



(14) 



where 



H^(r) =/!o(r) + U 
i2= J drl^,^'' 



(r) 



r(r)r + 2 



<At(r)A(r)> 



HJr) - m 



dt 



(15) 



This is effectively equivalent to the GPE in combination with a separate second-order 
correction |7|; the nonlinearity of the equations of motion prevents the expression of 
such a correction in closed form, however Note the nonlocal terms (off-diagonal forms 
of the normal and anomolous average), as well as the ultraviolet-divergent diagonal 
anomalous average, which must be renormalized lfT4l [TSl . The GGPE can also be 
generated from the approximate cubic Hamiltonian 

H2^H^^+Ao j dr'<A"l"(r')A(r')) 

= r dr[(/)^*ir) [2<A^(r)A(r))A(r) -i- A"'"(r)(A(r)2)] -i- H.c.} 



u 
u 



- ff drdr'[(l>^* 



(r)|0N(r)|2A(r) + H.c, 



(16) 



(r)|0''(r)|2 



X 



<A^(r')A(r))A(r') + A"'"(r')<A(r')A(r))] -i- H.c.}, 



^Note, however, Morgan's remark that Hartree-Fock-Bogohubov factorizations of third- and fourth-order 
terms neglect third-order corrections as large as fourth-order terms that are retained H4J . 
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in conjunction with the more complete form of the commutator [A(r), A^(r')] ~ <3(r, r')- 
(A' (r')A(r))/A^c- Both assume a Gaussian approximation |8|, and H^', M) are equiva- 
lent to H2 of Eq. A, with A^, U replaced by N^, tj. 

Substituting the GGPE back into the second-order equations of motion for A(r), 
A' (r) generated by Ht, eliminates all zeroth- and second-order terms: what remain 



(discarding higher-order terms) are the modified EDGE [Eq. ( 10 1], with U, N, A(r) 
replaced by U, Nc, A(r). Similarly, substituting the time-independent GGPE back into 
Hj, yields 

H^^H^'+Ao J flfr<Al(r)A(r)> - I J dr\<p''(r)f 

+ J fl'rA"'(r)[/io(r) + 2f7|0N(r)|2-^o]A(r) (17) 

where Hq \ is equivalent to Hq, with A^, U replaced by Nc, U, and we note that the first 
corrections to Hq - Hq' + '^o / dr{A ' (r) A(r)) appear at quartic order; other corrections 
Xo U - U are smaller than need be accounted for. H-^ and H2 are therefore effectively 
identical for a stationary state, and is literally identical to the stationary form of 
the second-order approximate Hamiltonian determined from an asymptotic expansion 
in powers of N^^^^ |8| (rather than A^"'^^ 16] |2l). In terms of quasiparticles bu, b\ 
[defined with respect to A(r), A' (r) as S^, P^, are defined with respect to A(r), A'^(r)] 
this becomes (assuming a thermal equilibrium state and bosonic quasiparticles): 



H, - IJ dv\<l>''{r)t + g 6, ^b[h - J c 

+ AoY,\dr [{{blh) + l)\vk{rf + {b[bu)[\vk{rf + 1]} . 



(18) 



The significance of deviations from exactly bosonic quasiparticle commutation rela- 
tions with increasing depletion is difficult to quantify in the abstract jH); we note that 
this does not seem to have been an issue in the finite-temperature calculations of Mor- 
gan, however [23 , 15. 24] . 

The A2 is in general complex; this constrains (p^ir) to unit norm even though 
the GGPE includes terms describing particle transfer between condensate and non- 
condensate. It may be advantageous to consider a condensate wavefunction with vary- 
ing norm: 0(r) - ^/NJN(p^(r). Noting that A^^ - N - Jt/r(A"'"(r)A(r)>, condensate 
number dynamics can be written (to quadratic order) as thdNddt - U J dr[(/)^*(ry^{A(r)^)- 
{AHrf}^^(rf] = (A2 - Al)Nc |8|. Hence, from md^{r)ldt = ^^NJN{m^^^ (r) I dt] + 
(mdNc/dt)^(r)/2Nc, 



/;z|<D(r) = [//g(r) - Ar] (5(r) + u^*(r)^-^^^ - ^ f dr'mr')f 

ut TV iSc J (19) 

X [<A"'"(r')A(r))0(r') + <t*(r')<A(r')A(r))] , 
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where 



HJr) =/i„(r) + U 



1 \,^._m2 , .<A^(r)A(r)) 



1- — |<l)(rr + 2 



Nc N 



(20) 

Here is real, and so, like A, simply describes an arbitrary phase. Carrying this 
formulation over to the modified EDGE, we note that, e.g., lJ\(p^{r)\^ - f/|4)(r)|^. 
Hence, the correct equations are obtained by setting 0^(r) — > <l)(r) in Eq. (iTOjl. Another 
advantage is that the condensate number can be tracked by A^c — J t/r|<l)(r)p7ratherthan 
the more involved procedure of continuously determining Nc - N - J dr{A^ (r)A(r)). 

2.5 Numerical Implementation 

Within the number conserving Bogoliubov (i.e., dynamically first-order) treatment, 
first determine an appropriate stationary solution ^'^(r) of the time-independent GPE, 
with its corresponding A. With these, construct and diagonalize X(r, r') to determine 
its eigenmodes. Choosing the quasiparticle operators to be time-independent |7|, the 
dynamics of A(r), A ' (r) are entirely determined by 

»fM. rrfr'£(r,r')("f,'M^ (21) 



(22) 



dt\ vk(r) j J \ vkir') 

propagating as many initial eigenmodes as deemed necessary. During the calculation, 
one can discard the projectors Q(r, r'), Q*{r, r'), regaining the standard EDGE. Apply- 
ing 

Q(r,r') 
Q*(r,r') 

to the final propagated modes then gives the correct result fT). One must propagate the 
time-dependent GPE (with A set - 0) in parallel, feeding the solution into the EDGE 
(and, ultimately, the projectors). Nonlocality of the modified EDGE is in this case not 
a practical issue. 

Determining a stationary GPE solution is also the starting point when applying the 
second-order treatment. To study dynamics starting from an ultracold sample with ini- 
tially negligible depletion, the initial setup should be equivalent, except that one must 
determine the normal and anomalous average terms to feed into the GGPE (note that 
the anomalous terms in particular are generally slow to converge, although semiclassi- 
cal approximations can significantly help 1 15 /24]). If the dynamics induce substantial 
depletion, the GGPE evolution will soon differ from that produced by the GPE. In the 
modified EDGE to which the GGPE is coupled, the projectors do not now separate out; 
hence, nonlocality due to the projectors and the off-diagonal normal and anomalous 
terms is a significant numerical issue. A finite temperature initial condition must be 
produced from an initial stationary GPE solution in an iterative, self-consistent manner 
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3 Validity Issues 



The treatments presented here are effectively perturbative expansions of increasing or- 
ders of fluctuation terms about a classical field (the condensate). As such, their validity 
relies upon {Na^)^''^ « 1 if T = 0, and {k^T I gN^){N^a^yi^ <s 1 if {k^TjgN^) » 1, 
where T is the temperature and k-^ is Boltzmann's constant f\A\. More stringently, 
dynamics propagated by the GPE coupled to the modified BDGE are only valid so 
long as the the non-condensate fraction [as determined by Jiir(A' (r)A(r))] remains 
insignificant, i.e., 1 - N^jN <K 1. 

The second-order Gardiner-Morgan treatment |8| is the basis of Morgan's analysis 
of excitations to finite-temperature BEC Il23l [TSi .24 J , to good agreement with exper- 
iment ll26l . Good agreement was also achieved also achieved by the ZNG treatment 
IZfj of Jackson and Zaremba |28|. We note this treatment (unlike Gardiner-Morgan 
|[8J) does not appear to explicitly account for the phonon character of low-energy 
states or the anomalous average and Beliaev processes ||23]| . which can be significant 
lfT4l |29l l30l [3n l32ll . In summary, ZNG seems more consciously oriented towards 
regimes of thermal equilibrium (and as such is numerically more tractable), whereas 
the approach presented here includes effects that are likely to be more significant at low 
initial temperatures or situations far from equilibrium. 

Another contrasting approach able to account for a significant non-condensate frac- 
tion is that based on a cumulant expansion developed by Kohler 1341 and coworkers. 
This has been very successful in describing the formation of Feshbach molecules, and 
also accounts for the dynamical loss of (non-molecule) atoms from the initial conden- 
sate fraction ll35l . An essential difference in this approach is the inclusion of a more 
physical scattering potential capable of supporting bound states (molecules), conse- 
quently also avoiding renormalization issues. As such, significant condensate depletion 
can occur due essentially to pure two-body, relatively high-energy scattering processes, 
without necessarily addressing loss due to the low-energy, many-body processes re- 
cently studied by Billam and Gardiner ||33i, and others l,2()l [3611371 [38., 2 1..39J. 

4 Applications 

Morgan's analysis of finite-temperature BEC excitations 1231 [TSl [24ll [see Fig. [TJ^)] 
involved applying a linear response treatment to the second-order equations, justified 
by the excitations being due to small perturbations; the equations take a rather involved 
appearance as a result, but are numerically more tractable than a full dynamical calcu- 
lation. A fully dynamical treatment, however, is be necessary to study non-perturbative 
dynamics [1201 1211 1391 where significant depletion from an initial very low tempera- 
ture Bose-Einstein condensate is expected, self consistently. This has recently been 
carried out by Billam and Gardiner within a simple quasi- Id 5-kicked-rotor-BEC con- 
figuration [see Fig. [TJb)]. The kicked rotor is well known system in the context of 
chaotic and quantum chaotic dynamics, and has had numerous atom-optical realiza- 
tions l40l[4n[42l[43l . as well as being an ideal dynamical test-system. An important 
conclusion is that a BEC appears to be impressively robust, in contrast to what might 
be expected from first-order treatments i[20l[37l[38l[21il . 
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Figure 1 : Implementations of the second-order self-consistent number-conserving ap- 
proach to Bose-Einstein condensate and non-condensate dynamics developed by Gar- 
diner and Morgan |8 1. (a) Figure adapted from |23 1 comparing experimental |26 1 (solid 
circles) and theoretical (open circles) excitation energies for excitation modes with ax- 
ial quantum numbers (i) w = and (ii) m = 2 as a function of reduced temperature 
T/T^. Here r'? is the BEC critical temperature for an ideal gas, and cOr is the radial 
trapping frequency, (b) Figure adapted from ||33]| . of sample dynamics of a 5-kicked- 
rotor BEC, a simplified spatially -periodic quasi- Id system with time-periodic driving 
potential V(6, t) - KCOs(ff} YjJ=[) ^(t ~ jTp), showing condensate and non-condensate 
n, fractions evolved by (i) first-order and (ii) second-order number-conserving de- 
scriptions. (b)(iii) plots of a coherence measure 

gi(0, 0') - {'¥'''' (6')'i'(6)}/N, and the fidelity of the condensate wave function to the GPE 
evolution F - \ ^ d0^*^^^(p{6)\^ for dynamics corresponding to (b)(ii). 
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